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^ : Abstract 
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t~^ ' We construct a driven sandpile slope model and study it by numerical sim- 

O" 

ulations in one dimension. The model is specified by a threshold slope ac, a 

■r-j- ' parameter a, governing the local current-slope relation (beyond threshold), 

Q\' 

and jin, the mean input current of sand. A nonequilibrium phase diagram is 
obtained in the a — j\^ plane. We find an infinity of phases, characterized 
rH ' by diff'erent mean slopes and separated by continuous or first-order bound- 

o' 

O ■ aries, some of which we obtain analytically. Extensions to two dimensions are 

\_^ • discussed. 
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The statistical mechanics of nonequihbrium steady states is a subject of growing general 
interest. Phase transitions between such states are by no means as well understood as 
their equilibrium counterparts. Some insight has been gained into this problem by the 
study of simple driven lattice models |]1],3. Sandpile models, which introduced the notion 
of self- organized criticality (SOC) [0] as a general explanation for the wide occurrence of 
power laws in nature, are a natural setting in which to study phase transitions far from 
equilibrium. Surprisingly, to our knowledge, this has not been attempted. We construct 
a driven sandpile model and show that it exhibits, in spite of its simplicity, a rich phase 
diagram, thus making it a good laboratory for the study of nonequilibrium phenomena. In 
particular, our model exhibits continuous transitions from pinned or threshold-dominated 
states to unpinned states; these are reminiscent of dynamical phase transitions in more 
complex systems such as sliding charge- density waves (CDW's) Q and pinned flux-lattices 

i- 

Earlier studies of sandpile models have concentrated on SOC, either in steadily flowing 
sandpiles [^J^ or at the angle of repose |Q. Much of the work has been on critical-height 
models, in which the update rule depends only on the height at each site; critical- slope models 
(CSM) [^ have been studied to a lesser extent. In this Letter, we present a comprehensive 
study of a simple CSM in which the current-slope relation for slopes exceeding a threshold 
is controlled by a parameter a. We monitor the steady states of our model as a function of 
a and the mean input current jin and find the rich nonequilibrium phase diagram (Fig. 1): 
It shows that there are many phases characterized by the average slope cxav of the sandpile. 
In the repose phase, a^.^ = cTc, the slope at the angle of repose. As jin is increased at fixed 
a, the repose phase undergoes a continuous transition (solid line) to an unpinned super- 
repose phase at which o"av — cTc rises continuously from zero with an exponent (3 c^ 0.5 (Fig. 
2). This continuous transition is followed by a series (which we argue is infinite) of similar 
continuous transitions. The repose phase lies between two first-order boundaries: one at low 
a to a pinned super-repose (cXav > o'c) region, the other at large a, to a sub-repose phase 
(cTav < cTc)- The lower one of these first-order lines meets the continuous line at jin = 0.5 at a 



multicritical point. The pinned super-repose region contains an infinity of phases, separated 
by first-order hneces parallel to the jin axis (Fig. 1); (Xav jumps at these boundaries. We 
also monitor height profiles, the equal-time height correlation function and the associated 
correlation length, and current autocorrelations and the associated correlation time. Both 
the correlation length and time (Fig. 3) diverge at the continuous transitions in Fig. 1. We 
show that our main results can be understood on the basis of a mean-field theory and a 
mode-softening argument. 

In our model, integers hi specify the heights of columns of sand at the sites i of a one- 
dimensional chain {1 < i < N). The stability of the column at a site is determined by a 
threshold condition which mimics the angle of repose for a real sandpile: When the height 
difference between a site i and its right neighbor {i + 1) (i.e., the local slope ai = hi — /ij+i) 
exceeds ac, some sand topples to the right neighbor, and hi is updated via 

hi~^ hi- ji , hi+i -^ hi+i + ji , (1) 

ji=X{axa,)e{a^-ac). (2) 

Af{x) is the integer nearest to x, B the step function, and a a real number. This part of 
our update rule conserves the number of particles locally and yields a local current which 
increases with the local slope, thus preventing an unbounded buildup of particles in the pile 
0. The parameter a controls the current-slope relation for slopes exceeding ac- We restrict 
our study to a > 0, since a < yields unphysical runaway behavior; the upper bound for 
a is chosen to limit the region we explore. The mean input current of sand particles ji^ is 
another control parameter. At each time step, we add m particles to a randomly chosen site 
with probability p, so jin = p y< m. We set m = 10 for specificity (our results do not depend 
on this choice) and cover the range 0.001 < p < 0.15, so 0.01 < ji^ < 1.5. This addition 
of particles violates local particle conservation; and, for such an addition rate, the mean 
input current and the noise amplitude per site vanish as A^ ^ oo. Particles are allowed to 
leave our system through the right, but not the left {i = 1), boundary: any particle that 



reaches the A^*'* site is removed immediately. Our boundary conditions and update rules 
(1-2) clearly pick a direction for the current (from left to right). 

We use initial conditions in which hi = ac{N — i) +6i, where 6i is an integer that assumes 
the values 0, ±1 with equal probability. We update all sites simultaneously and allow the 
system to evolve till it attains a steady state, i.e., when jjn = jout, the average output current 
(the average number of particles dropping out from the open boundary per unit time). In 
practice, we say that the steady state has been achieved when these two currents are within 
1% of each other. Once the steady state is obtained |l^, we accumulate data for 10^ — 10^ 



updates per site (UPS). Data for averages are stored after every 50 UPS. We also average 
our data over 10 — 50 different initial conditions. To minimize boundary effects we ignore a 
few sites (three or more if necessary) near each boundary while computing all averages. 

The nonequilibrium steady state of our model can be characterized by the mean slope 
o"av (the local slope ai averaged over i and many time steps), which is the order parameter 
for our model. We also monitor the mean current jav, the output current jout, the equal-time 
height correlation function Chhij) = {{[h\ ~ {^i)t\Wi+r " {^i+r)t\)i}v ^^^ ^^^ output current 
autocorrelation function Cjj{t) = ([j*^^ - {jout)t][fott ~ Oout)t])i, where (■ ■ ■)i and (■ ■ ■)i 
denote, respectively, averages over i and time t. 

In Fig. 2 the asymptotes (dashed lines) a"av = Oc and cTav = Jin/ (2a) indicate the behavior 
of CTav at very low (< 0.5) and very large (^ 0.5) jin, respectively. The full curve (for a = 0.1 
and A^ = 128 in Fig. 2) shows that the approach to these limits is nontrivial: As jin increases 
there are successive onsets, indicated by arrows. There might well be an infinity of such 
onsets (see below), but they become hard to resolve numerically at large jin. The inset shows 
a finite-size scaling plot at the first of these onsets. For jin ^ 0.5 , o"av = (Tc = 10, i.e., we 
have threshold-dominated behavior at low jin. The asymptotic behaviors can be understood 
via the mean-field theory presented below. 

The evolution equation for hj is 

ht'-K = -jl+jU + vU (3) 



where the noise rjj has a mean jirjN that accounts for the addition of particles. If we average 
over this noise, then, in the steady state, we get jj— jj-i = jm/N. If we impose the boundary 
condition ji=N = Jm, we get ji = ijin/N and thence a mean current jav = {{ji)i)t = ^Jm ~ 
jin/2 for A^ -^ cxD. We use these exact results to check our simulation. Our mean-field theory 
assumes that 

{{Ar{aa,)Q{a, - a,))J, ^ ((Ar(aa,))J, {{Q{a, - a,))J, . (4) 

For large jin, most (Tj > (Tc, so we further assume ((6((T.j — o'c))j)t — 1. If its argument is 
large, the discontinuities of A/'(acrj) are small relative to aai, so, on averaging Eq. (H), we 
make the approximation jav = {{■^{c(<^i))i)t = «crav, whence Uav ~ jin/(2a), the large-jin 
asymptote of Fig. 2. Given these approximations, Eq.(|^) yields a discrete diffusion equation 
for hj with a spatially uniform source jin/N. If we solve this with our boundary conditions 
we get parabolic height profiles for large jin, which we also find in our simulations [^ . 



The the current-slope relation Eq. (0), shows that (Fig. 4), for a uniform profile with 
slope (J, no current flows if cr < cXc- For a > cXc, the current grows with the slope in discrete 
steps, which reflect the Af function in Eq. (|^). As we increase a, the width of the step at 
ji = 1 shrinks till it becomes a single point at a = 0.125, the value at which the continuous 
transition at j-^ = 0.5 (Fig. 1) terminates. This can be understood as follows: If we turn 
Fig. 4 on its side, we see that the slope is pinned at 0"^ below some threshold value of ji. 
Beyond this threshold, the slope rises sharply before it saturates at another value of a. The 
onsets in Fig. 2 are just the sharp steps of Fig. 4 rounded by our spatiotemporal average. 

The vanishing of a step in Fig. 4 can also be linked with the termination of the continuous 
line via a mode-softening argument. Consider, e.g., the step at ji = 1 for which the 11 < 
ai < 14 (Fig. 4 with a = 0.1). All slopes o"j in this interval are equivalent in the sense that 
the sandpile dispenses the same amount of local current for all of them. Thus the restoring 
force, in response to a change of an onsite slope from a to o" + 1, vanishes and leads to 
divergent correlation lengths and relaxation times (Fig. 3). Clearly the infinity of steps in 
Fig. 4 imply an infinity of onsets in Fig. 2, though the large- jm onsets are hard to resolve 



numerically. 

The above arguments do not yield the value of jin at onset since our spatiotemporal 
average shifts the values of the thresholds in Fig. 4. The actual value of jin at the onset 
depends on the distribution of slopes in the interior of the pile. We have done some numerical 
and analytic calculations []TT| on a 'single-step model' 0, which justify the occurrence of 
the transition at jin = 0.5. 

For the first onset we calculate the mean slope near jin — 0.5 for a = 0.1 and A^ = 
32,64, 128, 256and512 and from a finite-size-scaling analysis (Fig. 2) obtain the exponents 
P = 0.5 and z/ = 4.0, where J = |jin — jc\/jc with jc = 0.5, and we use the scaling form 
Sa = N'^/^^F^J^N)^ with 5a = o"av — <^c- The exponent z/, which we obtain in this way from 
finite-size scaling, might not be the same as the actual correlation length exponent as has 
been pointed out in the context of sliding CDW's 0. 

In the low- jin regime, as a increases from to 0.05, (Tav decreases in steps of one. The 
values of a at which these steps occur yield the first-order boundaries of Fig. 1. All these 
first-order boundaries end at critical points in the range 1 ^ jin ^ 1.2. The overall envelope 
of the steps in o"av can be fit approximately to a form o"av ~ l/(2a). This behavior follows 
from the update rule (Eq. ^: A/'(ao"j) = for aai < |, so, if a < \/{2a^ (=0.05 for 
ac = 10), even a slope cij > ac may become stable; e.g., with ac = 10 and a = 0.04, 
N{aai) = for cTj = (Jc + 1 and cTc -|- 2, so the slope at the effective angle of repose turns out 
to be a(,Q = (Jc + 2 = 12. As a decreases, a^E also increases in discrete steps, leading to the 
first -order lines at low a. 



We fit the equal-time height correlation function Chh{f) to an exponential form ||TT 
which then yields a correlation length. Figure 3 shows this correlation length ^hh versus jin 
for a = 0.1 at different values of A^. Clearly ^^h diverges at jin — 0.5 and also at subsequent 
continuous transitions. We have not characterized the divergences by an exponent because 
^hh exceeds our system size somewhat before the transition (so we do not show data near 

Jin = O.5.). 

The output current autocorrelation function Cjj{t) does not fit an exponential form 



very well. However, we have extracted correlation times Tjj from the area under Cjj{t) 
(normalized so that Cjj{t = 0) = 1). Figure 3 shows plots of Tjj versus jhi for a = 0.1, 
which sharpen with increasing A^. This sharpening shows up clearly near jin = 0.5 and leads 
to an increasing trend in Tjj near the next onset (jm — 1.1). 

It is generally believed that nonequilibrium phase transitions cannot occur in stochastic, 
one-dimensional models with short-ranged interactions ||12|. We believe the transitions in 



our model occur because the noise amplitude per site vanishes as A^ ^ oo. This is why our 
one-dimensional model with short-range interactions exhibits phase transitions. We have 
checked that these phase transitions get rounded if there is a finite noise amplitude per 
site as A^ — > oo. We have also performed simulations on two-dimensional versions of our 
model. Our preliminary results indicate that, with low noise and a variety of boundary 
conditions, such two-dimensional models show the same transitions as the one-dimensional 
model discussed here. In the steady state, the two-dimensional system behaves like uncou- 
pled one-dimensional systems, so no transitions occur in the high-noise limit. The details of 
this study will be published elsewhere |]TT| . 

We have shown that our driven sandpile model displays a variety of steady states and 
many transitions between them. This richness, coupled with its simplicity, makes our model a 
promising one for the study of nonequilibrium phenomena in driven systems. As noted above, 
it displays transitions similar to those in other driven systems, e.g., our onset transitions 
(Fig. 2) are like unpinning transitions in sliding CDW's |^ or in pinned flux-lattice systems 
0. It would be interesting to study whether this similarity is merely superficial. There 
are some obvious ways in which the CDW models are different from ours: (1) They exhibit 
pinning because of quenched randomness but have no external noise; and (2) no current 
flows in their pinned states. The importance of these differences needs to be elucidated. In 
this general context it is interesting to study the zero-current limit of our model. We find 
that it does not show conventional SOC when the pile is allowed to relax, after each input 
of sand, to a completely quiescent state in which no further transfers are possible (a la Bak, 
Tang, and Wiesenfeld P|). The precise forms of the distribution of avalanche sizes, etc., will 



be reported elsewhere |TT[ . 
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FIGURES 
FIG. 1. The nonequilibrium phase diagram of our model in the a— Jin plane. Full (dashed) lines 

indicate continuous (first-order) transitions. At low a there is an infinity of first-order boundaries; 

only the first five are shown. The corners in the phase boundaries and multicritical points are 

schematic; our data are too noisy to resolve them. 

FIG. 2. A plot (full line) of cJav versus jin- Dashed lines are the asymptotes for large and small 
Jin- Arrows mark successive onsets (see text). The inset shows a finite-size scaling plot for baN'^l^ 
versus N r at a = 0.1 and jin > 0.5 for iV = 64(+), 128(o), 256(*) and 512(A) with /3 = 0.5 and 
i/ = 4. 

FIG. 3. Equal-time height correlation length ^/j/^ and output current autocorrelation time Tjj 
versus jin for N = 32(o), 64(+), 128(o), 256(*) and 512(A). For ^/^/j we do not plot points in regions 
where S,hh> N. 

FIG. 4. A plot of ji versus ai (Eq. (2)) for a = 0.1 and 0.125 and a^ = 10. 
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